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Polyharmonic  Splines  in  IRd: 
Tools  for  Fast  Evaluation 


R.  K.  Beatson,  J.  B.  Cherrie,  and  David  L.  Ragozin 


§1.  Introduction 

As  is  now  well  known,  hierarchical  and  fast  multipole-like  methods  can  greatly 
reduce  the  storage  and  operation  counts  for  fitting  and  evaluating  radial  basis 
functions.  In  particular,  for  spline  functions  of  the  form 

N 

s(x)  =p(x)  + ^2\i<l>(\x  -Xi\),  (1) 

i=l 


p a low  degree  polynomial,  the  cost  of  a single  extra  evaluation  can  be  reduced 
from  O(N)  to  0(1)  operations,  and  the  cost  of  a matrix-vector  product  (that 
is,  evaluation  at  all  centers)  can  be  decreased  from  0(N 2)  to  0(N). 

This  paper  outlines  some  of  the  mathematics  required  to  implement  meth- 
ods of  these  types  for  polyharmonic  splines  in  Rd,  d even,  that  is  for  splines  s 
corresponding  to  (j)  chosen  from  the  list 


j r2(f+i-d/2))  e = 0,...,d/ 2-2, 

( r2C+1-<i/2)  log(r),  C = d/2  — 1, ...  . 


(2) 


We  carry  out  most  of  our  work  in  the  general  Rd  setting  and  then  specialize 
to  d = 4.  We  refer  the  reader  to  our  more  detailed  work  [1]  which  contains 
all  the  details  of  this  special  case.  We  are  currently  working  on  developing  all 
the  details  for  the  general  Rd  case. 

A key  technique  in  our  development  is  the  exploitation  of  the  rotation 
group  invariance  of  radial  basis  functions.  This  means  that  we  exploit  the  fact 
that  any  kernel  k(x,  y)  = <f>( \x  — y\)  will  be  rotation  invariant  in  the  sense  that 


k(gx,gy)  = k(x,y),  for  all  orthogonal  g 6 0(d).  (3) 


Invariance  leads  to  many  crucial  simplifications  and  efficiencies  in  developing 
and  manipulating  the  polyharmonic  expansions  which  lie  at  the  heart  of  the 
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hierarchical  and  fast  multipole  expansions.  Related  development  of  general 
spherical  harmonic  expansions  based  on  these  techniques  can  be  found  in  [6, 
7,8]. 

We  will  not  detail  the  basic  framework  of  hierarchical  and  fast  multipole 
methods  within  which  this  mathematics  sits.  However,  we  do  recall  that 
an  essential  component  of  the  method  is  the  grouping  of  approximations  to 
summands  like  (1)  into  subsums,  which  are  approximations  to  the  influence 
of  that  part  of  (1)  associated  with  centers  in  a single  panel  or  cluster.  The 
key  steps  to  obtain  them  requires: 

• Finding  explicit  Taylor/Laurent  expansions 

For  each  x,  <p{\x  — £<!)  = E Pm(l,®<),  |z<|  < |*|,  (4) 

m>mo 

with  p m homogeneous  polynomials  of  degree  m in  z< . 

• Finding  an  efficient  separation  of  the  x,  z<  influence  in  , i.e.,  expanding 

(5) 

t 

for  some  good  choice  of  (basis)  functions  {/,m(z)}  and  {<?]"(£)}. 

These  expansion  and  separation  results  provide  the  approximations  to  sub- 
sums which  are  the  far  and  near  field  expansions.  Other  essential  components 
are  the  tools  to  manipulate  these  expansions,  namely  error  estimates,  unique- 
ness theorems,  and  translation  formulae.  In  this  paper  we  concentrate  on  the 
algebraic  tools  and  give  some  extensive  general  results  on  (4)  (Theorem  6) 
and  (5)  ((20)  in  Section  3)  and  for  IR4  we  give  the  appropriate  far  and  near 
field  expansions  for  the  <f>i  (Theorems  7 and  8),  and  a brief  indication  of  the 
dual  basis  leading  to  (20).  Analogous  results  for  polyharmonic  splines  in  IR2 
appear  in  [3].  The  reader  unfamiliar  with  the  framework  of  the  fast  multipole 
method  may  wish  to  refer  to  the  original  paper  of  Greengard  and  Rokhlin  [4], 
or  to  the  introductory  short  course  [2]. 


§2.  Polyharmonic  Functions  and  Homogeneous  Polynomials 

First  we  record  some  detailed  facts  relating  to  the  Laplacian  A,  and  its  actions 
on  special  homogeneous  functions  and  the  logarithm  of  the  distance.  In  par- 
ticular, we  show  why  our  basic  functions  (j)i  in  (2)  are  polyharmonic  or  more 
specifically  (£  + l)-harmonic  in  the  sense  that  A e+1<f>e  = 0. 

Lemma  1.  Let  \ ■ \ be  the  2-norm  on  3Rd,  d even. 

i)  If  f : ]Rd  \ {0}  — > IR  is  a non-trivial  harmonic  function  that  is  homoge- 
neous of  integral  degree  m,  then 


A(|  • 1 2<7)  = 2*(d+2£  + 2m-2)|  - l2^/- 


(6) 
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Hence  \ ■ \2tf  is  polyharmonic  of  exact  order 

( t + 1,  for  £ >0 , m > — | or  m < 1 - l—  | , 

\ l + m + | , for  £ < 0,m  > 1 — £ — 


ii) 


In  particular 


( {£  + l)-harmonic  for  £ > 0, 

( (£  + d/2)-harmonic  for  —d/2  < £ < 0. 


| • \2t  log  | ■ | is  (£  + d/2) -harmonic  for  £>0.  More  generally, 


A|  - 1"  log  | • | = | ■ \^t~1\2£{2£  + d - 2)  log  | • | + (4£  + d - 2)). 


(7) 


Proof:  For  the  first  part  of  (i),  just  apply  the  product  rule  for  the  Laplacian, 


A (fg)  = (A f)g  + 2 (V/)  • (V9)  + f(Ag), 


and  the  Euler  relation  for  a function  / that  is  homogeneous  of  degree  m, 

x ■ (V/)(x)  = m/( x). 

Then  observe  how  many  applications  of  A are  required  to  reduce  one  of  the 
multipliers  to  0.  Specializing  to  the  case  / = 1 yields  the  last  result  of  (i). 

The  first  part  of  (ii)  follows  from  (7)  in  combination  with  (i)  and  its  proof. 
(7)  follows  from  (i),  the  product  rule  for  the  Laplacian,  and  the  computation 
of  V log  | • | and  A log  | • | . □ 

From  the  detailed  eigenvalue-like  information  on  the  Laplacian  map  in 
(6),  we  can  get  a decomposition  theorem  for  II„,  the  homogeneous  polynomials 
of  degree  n,  in  terms  of  the  spherical  harmonics  of  degree  n: 

B„  = {p  € n„  : Ap  = 0}  = ker(A)  n IIn. 

This  is  useful  in  understanding  the  structure  of  the  homogeneous  polynomial 
terms  in  any  Taylor /Laurent  type  series  expansions  for  the  4>(\x  — Xi\).  In 
view  of  Lemma  1 the  decomposition  splits  IIn  into  its  harmonic,  biharmonic, 
triharmonic,  etc.,  parts. 

Lemma  2.  IIn  = | • |2<Hn_2£.  In  particular,  H„  n | • |2nn_2  = {0}. 

Proof:  Note  that  for  n = 0, 1,  IIn  = B„,  so  the  base  for  an  inductive  proof 
is  true.  Assume  the  decomposition  for  n — 2,  some  n > 2,  so  for  each  p G Hn, 

L(n-2)/2J 

A p=  ^ | • 1 2lhn-2-2i,  some  hi  G H n_2_2*. 

e=o 

Then  by  (6),  if 

l(n-2)/2j 

hn=P~  9(£+1)~1(d  + 27l-2(£  + 2)r1|  ' \2(-t+1)hn-2(t+l), 

e=o 

then  Ahn  = 0.  So  hn  € Hn  and  the  decomposition  of  IIn  is  proved  by  induc- 
tion. □ 
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Some  additional  consequences  of  (6)  come  when  we  study  what  happens 
for  negative  £.  Here  we  have  noted  that  | • (2-2,n~df  js  harmonic  whenever  / 
is  m-homogeneous  and  harmonic.  But  bringing  the  factor  of  (|  • |~2)m  inside 
the  m-homogeneous  function  / shows  | • |2-d/(7l  • |2)  is  harmonic  for  any 
homogeneous  harmonic  function  /.  In  fact  this  construction  is  independent 
of  the  homogeneity  order  of  /.  In  general  the  Kelvin  transform,  defined  by 
Kf(x ) = |x|2_d/(3;/|a;|2),  which  arises  from  inversion  in  the  sphere  followed 
by  multiplication  by  |rc|2— d,  maps  harmonic  functions  to  harmonic  functions. 
On  nn  Kf  = | . 1 2-d-2nf'  Associated  with  the  Kelvin  transform  are  the  spaces 
of  negative  degree  2 — d — m 

oo  oo 

0 I • \2eKMm+2e  = 0 I . \2~d~2m~‘2eMm+2i,  (8) 

«=0  t= 0 

which  are  useful  for  analysing  Laurent  (far  field)  series.  An  application  of 
Lemma  1 shows  that  Equation  (8)  displays  the  space  under  consideration 
split  into  its  harmonic,  biharmonic,  triharmonic,  etc.,  parts. 

§3.  Rotation  Invariance  and  Simplified  Taylor  Expansions 

The  decompositions  of  polynomial  spaces  in  the  previous  section  already  sim- 
plify the  Taylor/Laurent  type  expansions  (4)  we  need  to  determine.  To  make 
further  progress,  we  want  to  exploit  the  rotation  invariance  of  <j>( \x  — aq|). 
When  we  come  to  combine  subsums  in  (1),  we  will  want  to  fix  x and  concen- 
trate on  rotations  (orthogonal  matrices)  which  fix  x.  When  we  are  given  a pole 
p € lRrf,  we  let  Gp  = {g  : g € 0(d),  gp  = p}  denote  the  rotations  about  the  ray 
through  p.  So  the  function  f£(xK)  = <p(\x  — x<|)  satisfies  f^(gxK)  = /^(k<), 
for  all  g € Gx.  We  refer  to  any  function  / which  is  unchanged  by  rotations  in 
Gp  as  a p-zonal  function.  In  particular  we  have  the  p-zonal  harmonics 

= {h  € Hn  : h(gy)  = h(y)  for  all  g G Gp },  (9) 

and  the  p-zonal  homogeneous  polynomials  n£.  Now  the  Taylor/Laurent  ex- 
pansion of  /£  as  in  (4),  will  have  p^(x,gx..)  = p^(gx,5x<)  = p^(x,x<),  for 
g £ Gx  since  the  homogeneous  terms  must  remain  unchanged  under  rotations 
(see  Theorem  6).  Thus  these  terms  will  be  x-zonal  polynomials  as  a function 
of  xK.  What  is  the  general  stucture  of  nj  and  H^? 

Theorem  3.  Fix  a pole  x € ]Rli\{0}  . Let  Xo(')  = E Xi(')  = 2(x,  •).  Then 
there  exist  a unique  set  of  constants  am,  m > 1,  such  that,  the  inductively- 
defined  sequence  of  homogeneous  polynomials, 

Xm  + l =XTXm-am  + lM2H2Xm  — 1>  m>  1,  (10) 

consists  of  harmonic  functions.  Moreover, 
i)  Xm  ‘s  an  m-homogeneous  x-zonal  harmonic  function,  which  is  rotation 
invariant  in  the  sense  that  Xm(d')  = Xm(')  f°r  9 e 0(d). 
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Xm(x<)  = 0 ^ x,x<  G JRd,  (11) 

and  the  Xm(x<)  are  also  homogeneous  and  harmonic  as  functions  of  x. 

iii)  If  x 0,  {|  • |2/Xm-2 <(■).  i ~ 0,  • • • , Lm/2J}>  form  a basis  for  nm-  Id 

particular,  Xj  is  the  unique  (up  to  a scalar  multiple)  element  of  11) . 

iv)  For  m,  n nonnegative  integers,  the  kernels 

W,n  (*,*<)=  I* \2t+m~K\x<\2e+n-KxxK-2i(x<),  « = min  (m,n), 

m e n mod  2,  i = 0, . . . , [k/2J,  ^ 

form  a basis  for  the  space  of  all  rotation  invariant  polynomial  kernels, 
Pm,n{x,  x<)  which  are  homogeneous  of  degree  m in  x,  and  degree  n in  xK . 

Proof:  The  proof  of  the  existence  of  am  and  (i),  (ii),  and  (iii)  is  by  induction. 
For  m = 0 (i)  and  (iii)  are  trivially  true.  Let  0 ^ h G H*.  Then  h has  the 
form  h(-)  = c(p,-).  Since  h is  x-zonal  c(p,  ■}  = c(p,g-)  = c{g~1p,-),  for  all 
g € Gp.  This  implies  p has  the  same  direction  as  x.  Hence  h is  a multiple  of 
Xi  and  (i)  and  (iii)  follow  for  rn  ==  1. 

Now  induction  shows  that  (10)  defines  m-homogeneous  x-zonal  polyno- 
mials which  are  m-homogeneous  in  x,  for  any  choice  of  am+i . Also  they  are 
rotation  invariant.  To  complete  the  inductive  step  for  (i)  with  a fixed  x we 
need  only  show  that  there  is  a unique  am+i  that  makes  Xm+i(')  harmonic. 

From  the  homogeneity  in  x,  we  may  assume  |x|  = 1.  Since  X\Xm  is 
a homogenous  polynomial  of  degree  m + 1,  Lemma  2 asserts  that  there  ex- 
ist unique  homogenous  harmonic  polynomials  gm+i_2£  such  that  XiXm  = 
2feo+1^  I ' Since  Vxf  = 2xT , the  product  rule  for  the  Lapla- 

cian  and  the  inductive  assumption  that  Xm  is  harmonic  show  that 

A(XlXm)  = 4xT  • VXm  = 4^iXmi  1*1  = 1,  (13) 

where  dx  denotes  the  directional  derivative  in  the  (fixed)  direction  x.  Since 
AdxXm  = dxAXm  = 0,  it  follows  that  x?Xm  is  bi-harmonic  and 

A(XiXm)  = A(!  ■ I29m-i)  = 2 (d  + 2(m  - l))gm_i. 

Since  A maps  x-zonal  functions  to  x-zonal  functions  it  follows  that  qm-i  G 
MTrn_i  and  therefore  by  part  (iii)  of  the  inductive  hypothesis  qm-i  is  a multiple 
of  Xm— i-  Thus  the  existence  and  uniqueness  of  am+i  making  Xm+i  harmonic 
is  proved. 

We  now  turn  to  the  inductive  step  in  the  proof  of  (ii).  Using  the  rotation 
invariance  part  of  the  inductive  hypothesis,  Xm+i(ff-1')  is 

Xl  GrMXmGT1-)  - aro+lM2|S_1  ■ |2Xm-l(ff-1-)  = 

xTxSr-ora+iM2|-|2xSLi. 


(14) 
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Since  rotations  and  A commute,  the  left-hand  side  of  the  above  is  harmonic. 
Thus  the  right-hand  side  of  (14)  equals  Xm+ 1 ar>d  Um+i  is  independent  of  g. 
Using  homogeneity  it  is  also  independent  of  |x|.  Hence  am+\  is  independent  of 
x.  The  symmetry  in  x,  x<  of  (10)  then  implies  (11),  and  hence  the  homogeneity 
and  harmonicity  of  Xm+i(x<)  as  a function  of  x. 

We  now  turn  to  the  inductive  step  in  the  proof  of  (iii).  Since  J\Trn  = II m ^ 
and  xT  — «kXk  we  may  assume  |x|  = 1.  It  suffices  to  show  that  dimn^+1  < 
|m/2j  + 1 since,  by  Lemma  1,  {|  • |2fXm+i-2t>  ^ — 0, [m/ 2J}  is  an  inde- 
pendent set  in  n*+1.  Since  we  can  rotate  e\  to  x by  some  orthogonal  map, 
which  will  isomorphically  map  IT)j1+1  to  n^+1,  we  prove  our  dimensionality 
statement  for  x — ej.  To  analyze  the  value  f(y)  of  any  ei-zonal  function 
/,  we  choose  orthogonal  g±  € Gei  which  transform  y into  the  coordinate 
plane  spanned  by  the  first  two  basis  vectors.  The  two  possible  values  for  the 
transformed  y are 

9±V  = (y,ei)ei  ± \/|2/|2  - (l/.e i)2e2. 

Then  f(y)  = f(g±y)  = f{yi,±\J\y\2  - y2,0,...,0).  In  particular,  / must  be 
even  in  its  second  variable.  If  / € n|)J . j, 

L(m  + 1)/2J  . s 21 

f(y)  = cty?+1~2' l(  \J\y\2-y2i)  I for  some  Q. 

Hence  the  functions  y™+x~2l(\y\2  - yj)e  span  n^|+1. 

For  (iv)  we  just  note  that  pmin(x,  •)  = \x\mpm:n{x/\x\,  ■)  € n by  the 
homogeneity  assumptions.  Thus  the  basis  facts  from  (iii)  imply  there  are 
functions  bn^{x/\x\)  with 

ln/2j 

Pm,n(x>-)=  Yl  Kt{xl\A)\x\m\-\nX*J-2t- 
e=o 

The  rotation  invariance  of  pm  n and  of  the  terms  |x|m|  • |n_2fX2^X'  implies 
bn,l{gxl\gx\)  = 6n^(z/|x|)  for  all  rotations  g.  Rotating  x/\x\  to  e\  shows 
bn,e{x/\x\)  = b„}e(ei),  i.e.,  the  bU]t  are  constants.  Moreover,  the  homogeneity 
in  x of  x]  shows 

L„/2j 

Pm,n(*,  •)=  E Kl\AU-{n-m)\  ■ \2tXXn-2  f (15) 

1=0 

Since  the  left  hand  side  is  a polynomial  of  degree  m in  x , and  the  | ■ |2*Xn-2t 
are  independent,  each  \x\2l~^n~m')Xn-2i  associated  with  a nonzero  coefficient 
must  be  a polynomial  of  degree  m in  x.  Hence,  n — m = 2 j must  be  even. 
Also,  applying  the  second  part  of  Lemma  2,  1(1  — j)  = 2 t—(n—m)  >0.  If 
m > n the  proof  of  (iv)  is  done.  If  m < n,  then  reindexing  the  sum  in  terms 
of  (i  — j)  yields  (12).  □ 

The  following  result  is  known  [5],  but  is  included  for  the  sake  of  com- 
pleteness. 
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Theorem  4.  Define  a rotation  invariant  inner  product  (pairing)  for  functions 
on  the  unit  ball  B C Rd  by 

[f,h]=cv  f f(y)h(y)dy,  Cy1  = voJ{|y|  < 1}.  (16) 

h t»l<i} 

i)  If  f,  h are  homogeneous  of  degrees  m,n,  respectively  with  fh  E L1(B), 
then  [/,  h)  = 0 if  and  only  if  /f|s|=1j  f(y)h{y)dA  = 0,  i.e.,  the  integral  of 
fh  over  the  sphere  Sd_1  is  zero. 

ii)  If  m f-  n,  then  | • |2lHn  and  | • |2-*lHm  are  orthogonal  with  respect  to  this 
inner  product(pairing),  provided  m + n + 2(i  + j)  > —d. 

Proof:  For  (i)  just  introduce  polar  coordinates  (r,  y)  E [0,1]  x 5d_1.  Then, 
by  homogeneity  and  scaling  properties  of  the  area  of  {x  : \x\  = r}, 


[f,h]=cy  f rm+n+d  1 f f(y)h(y)dAdr,  m + n + d>  0, 

Jr- 0 JS d-* 

and  the  result  follows.  By  (i)  it  suffices  to  prove  (ii)  when  i = j = 0,  since 
integrals  of  a product  | • 1 2(l+Afh  on  Sd_1  do  not  depend  on  i,j.  Let  / € 
Hm  and  h € Hn.  Then  by  Green’s  Theorem  and  the  Euler  relation  for 
homogeneous  functions 


0 = [ (f{y)Ah(y)-h(y)Af(y)\dy 

= / (f{y)'Vh(y)-h(y)Vf{y)\-ndA=[  (n  - m)f(y)h(y)dA. 

J{\v\=i}  \ J J{\y\=A 


Thus  /[|9|=1}  f(y)h(y)dA  — 0,  and  the  analogous  relation  holds  for  integration 
over  the  ball  by  part  (i).  □ 

An  application  of  the  above  gives 

Lemma  5.  The  constants  am,  m > 2,  in  the  3-term  recurrence  (10)  defining 
the  x-zonal  harmonics  \m  °f  Theorem  3 are  positive. 

Proof:  By  Theorem  4 and  (10), 

0 = [Xm+l.Xm-ll  = [xfXm,Xm-l]  -«m+l[|  ' |2Xm-l.Xm-l] 

= [Xm> XlXm— 13  am+l[i  ' PXm-l.Xm-l] 

= [Xmi  Xm  + «rn|  ' |2Xm- 2]  “ «m+l[|  • PXm-l-Xm-l] 

= [Xm.Xm]  -Om+l[|  ' fXm-l.  Xm-l]- 


Hence,  am+1  = [x™,  Xml/0  ' |2Xm-i, Xm-il  > 0.  □ 

Now  part  (iii)  of  Theorem  3 leads  quite  directly  to  the  structure  of  near 
and  far  field  expansions  of  general  rotation  invariant  kernels  iplxjX^.  The 


54 


R.  K.  Beatson,  J.  B.  Cherrie,  and  D.  L.  Ragozin 


heuristic  that  a far  field  expansion  of  with  respect  to  x can  be  found 

from  a Taylor  expansion  with  respect  to  x<,  has  been  known  to  us  for  some 
while.  Theorem  6 below  gives  a proof  that  the  underlying  idea  is  correct 
inimportant  special  cases.  In  fact,  we  have  the  following  result  for  such  ip 
which  are  jointly  homogeneous  (ip(ax,  ax<)  — a2nip(x,  x<))  for  some  even 
integral  power  and  are  analytic  about  xK  = 0,  such  as 

V’nOc.Zc)  = \x  — x<|2n(log(|a;  — x<|2)  — log(|a;|2)).  (17) 


Theorem  6.  Let  ip(x,  xK),  x,xK  G lRri  be  rotation  invariant,  jointly  homo- 
geneous of  degree  2 n and  analytic  in  x,x<!  for  |x<  | < \x\.  Then  there  exist 
constants  c7pn  t such  that  the  Taylor  expansion  of  ip  about  x<  = 0 has  the  form 


oo  |^m/2j 

#*,*<)=£  E 

m= 0 1=0 


- 2£ 


(*<) 


oo  [m/2J 

= E E c",<|x|2(n+<-m)|x<|2V^<_2<(x). 

m= 0 i=0 


(18) 

(19) 


When  ip  is  (k  + 1) -harmonic  in  x<:  the  upper  limit  on  i in  (18)  or  (19)  is 
min{A;,  [m/2j }.  If  ip(x,  0)  = 0 then  the  lower  limit  on  m in  (18)  or  (19)  is  1. 

Proof:  The  terms  p))j(a:,a:<)  in  (4),  the  Taylor  series  of  ip(x,  x<)  with  respect 
to  x<t  are  degree  m homogeneous  polynomials  in  x<.  When  any  Taylor  se- 
ries is  grouped  by  homogeneity  with  respect  to  x<t  each  group  is  uniquely 
determined.  Since  only  the  term  p)0(pa:,px<)  in  the  series  for  ip^gx^x^  has 
homogeneity  m in  i<,  the  rotation  invariance  implies  that  pPP  is  also  rota- 
tion invariant.  Similarly  the  joint  homogeneity  of  ip  yields  pjP(ax,  ox<)  = 
a2nPrn.{x,  £<)■  Since  for  any  x,x<  there  is  a rotation  g (or  reflection  if  d = 2) 
which  interchanges  the  rays  through  x and  x<,  i.e.,  g(x/\x\)  = (i</|x<|)  and 
g(x</\x<\)  — (x/\x\),  it  follows  that 


M 2(m  n)Pm{ x,xK)  = \x\2(-m  n)pm(|x|x</|x<|,|i<|x/|x|) 


1 2m 


, Pm 


= M2(m  n^Pm(x<,  x)- 


Since  the  final  right  side  in  this  string  of  equalities  is  an  m-homogeneous 
polynomial  in  x,  we  see  that  the  terms  in  (4)  have  the  form  |x|2(n_m^pm(a;,x<) 
with  pm  a rotation  invariant  m-homogeneous  polynomial  in  each  of  x,xK. 
Hence  (18)  follows  by  Theorem  3.(iii).  □ 

The  separation  properties  in  (5)  can  now  be  achieved  by  further  use  of 
rotation  invariance.  Each  of  the  subspaces  | ■ 1 2lMj,  j + 2£  = n,  which  occur 
in  the  decomposition  of  nn  is  rotation  invariant.  Hence  it  has  a (unique) 
rotation  invariant  reproducing  kernel 


dim  Hj 

k(x,y)  = \x\2t\y\2e  E fi(x)fi(v ). 


*=0 


(20) 


Fast  Evaluation  of  Polyharmonic  Splines 


55 


where  {|  • 1 2tfi}  and  {|  • \2ifi}  are  any  bases  for  this  subspace  which  are  bi- 
orthogonally  dual  with  respect  to  some  rotation  invariant  inner  product,  e.g. 
the  inner  product  (16)  (see  [8].)  But  by  (12)  in  Theorem  3.(iii),  since  k(x,y ) 
is  (exactly)  ( l + l)-harmonic  as  a function  of  y.  and  is  homogeneous  of  degree 
21  + j in  both  x,  y, 

k(x,y)  = cjlt\x\2t\y\2eXj(y),  (21) 

for  some  Cjy  > 0.  Equating  (20)  and  (21)  provides  separation  of  the  influence 
of  x,  x<  in  the  expression  for  Xj(x<)-  A consequence  is  separation  of  x,  x< 
in  the  far  and  near  field  expansions  given  by  Theorem  6,  thus  allowing  the 
combination  of  the  expansions  for  several  centers  x<  = Xi,i  = 1, . . . , N,  into 
one  expansion  about  0. 


§4.  Expansions  in  R4 

In  this  section  we  use  the  results  of  Section  3 in  the  R4  case  to  outline  the 
explicit  expansion  formulae  for  the  4n  in  (2)  for  d = 4.  We  start  with  the  far 
field  expansion  of  the  potential  function  \x  — x<]-2. 

Theorem  7.  For  x,x<  € R4  with  |x<|  < |as|, 

OO 

\x-x<\~2  = ]T  |x|-2(*n+l)CmXx<(a,))  Cm  _ I (22) 

m= 0 

Proof:  Since  \x  — x<  |-2  is  harmonic  in  R4,  an  expansion  of  this  form  holds 
for  some  constants  cm  by  Theorem  6.  Using  Xm  (x)  = Xm(x<) > multiplication 
by  \x  - a;<|2  = \x\2  + |x<|2  - Xi(*<)  yields 

OO 

1=  E (cm-Cm_!)|xr2mX^<) 

771=0 

+ (Cm-2  - cro_iam)|x|_2ira-1)|x<|2x^_2(a:<), 

when  the  recurrence  (10)  is  used  and  the  geometrically  convergent  expansion 
is  rearranged  to  group  terms  of  common  homogeneity  in  xK.  Then  equating 
coefficients  using  (iii)  of  Theorem  3 shows  co  = 1,  cm  = cm_i  and  cm_2  = 
Cm-iOm • These  must  be  consistent  so  am  = cm  = 1 for  all  m.  □ 

We  now  outline  the  expansion  of  ipn(x,  x<)  from  (17).  This  gives  us  the 
bulk  of  the  far  field  expansion  for  0n+2. 

Theorem  8. 

oo  min{n+l,Lm/2j} 

^„(x,x<)=£  £ c^lxl^-^lx.l^l.^x),  (23) 

771=1  £=0 

where  the  non-zero  coefficients  ( are  given  by  the  formulae  , 

and  the  recurrence  c”+/  = c"  , — c”  , , — c"  , , , + c"  , „ , . 

771,  t m,l  771—1,7  771— l,c— 1 771  — Z, 7—1 

Proof:  The  form  of  all  the  expansions  follow  from  (18),  since  ipn(x,  0)  = 0. 
The  explicit  determination  of  the  c^,  t,  the  n = 0 case,  is  done  in  Lemma  4.4 


56 


R.  K.  Beatson,  J.  B.  Cherrie,  and  D.  L.  Ragozin 


of  [1].  The  recurrence  for  the  t follows  as  in  Theorem  7 from  (10)  with 
am  = 1 upon  multiplication  of  the  ipn  case  by  \x  - x<|2.  The  details  are  in 
Lemma  4.6  of  [1],  □ 

The  explicit  construction  of  bases  for  Hy  (and  dual  bases)  which  are 
needed  for  the  separation  results  can  also  be  significantly  simplified  by  use  of 
the  rotation  invariance  perspective.  A detailed  development  in  ]R4  is  in  our 
previously  cited  work. 
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